home *** CD-ROM | disk | FTP | other *** search
/ Libris Britannia 4 / science library(b).zip / science library(b) / PROGRAMM / PASCAL / 0514.ZIP / CRAYZ15.ARC / VDIDR.PAS < prev    next >
Pascal/Delphi Source File  |  1986-08-01  |  6KB  |  153 lines
  1. { Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, NY 13730 }
  2.  
  3. program main ;
  4.  
  5.      { Program:                                            }
  6.      {                                                     }
  7.      {    LINPACK SGEDI Test Driver.                       }
  8.      {                                                     }
  9.      { Version:                 Date:                      }
  10.      {                                                     }
  11.      {    1.5/TURBO Pascal 3.0     08/02/86                }
  12.      {                                                     }
  13.      { Description:                                        }
  14.      {                                                     }
  15.      {    Uses LINPACK SGECO to compute RCond for matrix   }
  16.      {    A and LINPACK SGEDI to compute determinant and   }
  17.      {    inverse of matrix A where A is set up as a       }
  18.      {    Hilbert matrix of specified order using SYSGEN.  }
  19.      {    Uses SGEFA and SGEDI to invert inverse to com-   }
  20.      {    with original.                                   }
  21.      {                                                     }
  22.      { Author:                                             }
  23.      {                                                     }
  24.      {    Adam Fritz                                       }
  25.      {    133 Main Street                                  }
  26.      {    Afton, New York 13730                            }
  27.      {                                                     }
  28.  
  29. {-I DizZ.con                    CONSTANT Declarations      }
  30. {-I Dizz.typ                    TYPE Declarations          }
  31. {-I Dizz.var                    VARIABLE Declarations      }
  32. {$I CrayZ.con                   CONSTANT Declarations      }
  33. {$I Crayz.typ                   TYPE Declarations          }
  34. {$I Crayz.var                   VARIABLE Declarations      }
  35.  
  36.      i, j           : integer ;
  37.      aaID           : vARRAY ;
  38.      pID            : vARRAY ;
  39.  
  40. {-I DizZ.pas                    Vector I/O Routines        }
  41. {$I DrivZ.pas                   Vector I/O Routines        }
  42.  
  43. {-I CGen.pas                    Test System Generator      }
  44. {$I Hilgen.pas                  Test System Generator      }
  45. {-I VectScal.p                  MathPak (C) Routine Package }
  46. {-I SkipVS.p                    MathPak (C) Routine Package }
  47. {-I mpBLAS.pas                  MathPak (C) BLAS           }
  48. {$I BLAS.pas                    Basic Linear Algebra       }
  49. {$I vSGETP.pas                  Virtual Array Transpose    }
  50. {$I vSGEFA.pas                  LINPACK Factor             }
  51. {$I vSGEDI.pas                  LINPACK Determinant and    }
  52. {                                  Inverse                 }
  53. {$I vOUT.pas                    Virtual Array Output       }
  54. {$I OUT.pas                     SICE Output Routine        }
  55.  
  56. begin
  57.                                 { Initialize }
  58.    writeln('LINPACK SGEDI Test Program, CrayZ Version 1.5.') ;
  59.    writeln ;
  60.                                 { Get Order }
  61.    n := 0 ;
  62.    while (n < 1) or (n > lda) do begin
  63.       write('Order: ') ;
  64.       readln(n)
  65.    end ;
  66.                                 { Allocate Original Matrix }
  67.    vCreate (aID,'aMATRIX.$$$',n) ;
  68.                                 { Get Print Code }
  69.    write('Print Code: ') ;
  70.    readln (PrintCode) ;
  71.                                 { Generate Test System }
  72.    SYSGEN (aID, lda, n, b) ;
  73.    if PrintCode > 0 then begin
  74.       writeln ;
  75.       writeln('Original System (by column):') ;
  76.       writeln ;
  77.       vOUT (aID, n) ;
  78.       OUT (b[1], lda, n, 1)
  79.    end ;
  80.                                 { Allocate Copy Matrix   }
  81.    vCreate (aaID,'aaMATRIX.$$$',n) ;
  82.                                 { Fill Transpose Matrix.  }
  83.    for i := 1 to n do
  84.       Aj[i] := 0.0 ;
  85.    for j := 1 to n do
  86.       VectorWrite (aaID,n,1,j,n,Aj) ;
  87.                                 { Transpose Original }
  88.    vSGETP (aID, aaID, n) ;
  89.    if PrintCode > 0 then begin
  90.       writeln ;
  91.       writeln ('Transpose (by column):') ;
  92.       writeln ;
  93.       vOUT (aaID, n)
  94.    end ;
  95.                                 { Factor Matrix }
  96.    vSGEFA (aID, lda, n, IPvt, InfoCode) ;
  97.    if PrintCode > 0 then begin
  98.       writeln ;
  99.       writeln('Factored Matrix (by column):') ;
  100.       writeln ;
  101.       vOUT (aID, n)
  102.    end ;
  103.                                 { Compute Determinant and Inverse }
  104.    JobCode := 11 ;
  105.    vSGEDI (aID, lda, n, IPvt, Det, Work, JobCode) ;
  106.    if PrintCode > 0 then begin
  107.       writeln ;
  108.       writeln('Inverse (by column):') ;
  109.       writeln ;
  110.       vOUT (aID, n) ;
  111.       writeln('Determinant: ',Det[1]:14:8,'E',Det[2]:3:0)
  112.    end ;
  113.                                 { Allocate Product Matrix }
  114.    vCreate (pID,'pMATRIX.$$$',n) ;
  115.                                 { Original Times Inverse }
  116.    for j := 1 to n do begin
  117.       iAj := VectorRead (aaID,n,1,j,n,Aj) ;
  118.       for i := 1 to n do begin
  119.          iAi := VectorRead (aID,n,1,i,n,Ai) ;
  120.          Ak[i] := SDOT (n, Aj[1], 1, Ai[1], 1)
  121.       end ;
  122.       VectorWrite (pID,n,1,j,n,Ak)
  123.    end ;
  124.    if PrintCode > 0 then begin
  125.       writeln ;
  126.       writeln('Original times Inverse (by column):') ;
  127.       writeln ;
  128.       vOUT (pID, n)
  129.    end ;
  130.                                 { Try to Restore Original }
  131.    vSGEFA (aID, lda, n, IPvt, InfoCode) ;
  132.    if InfoCode = 0 then begin
  133.       vSGEDI (aID, lda, n, IPvt, Det, Work, JobCode) ;
  134.       if PrintCode > 0 then begin
  135.          writeln ;
  136.          writeln('Inverse of Inverse (by column):') ;
  137.          writeln ;
  138.          vOUT (aID, n)
  139.       end
  140.    end
  141.    else
  142.       writeln('Error: Inverse is Singular.') ;
  143.                                 { Close Virtual Arrays }
  144.    vClose (pID) ;
  145.    vClose (aaID) ;
  146.    vClose (aID) ;
  147.                                 { Done }
  148.    writeln('End of Test.')
  149.  
  150. end.
  151.  
  152. { Copyright (C) 1986 Adam Fritz, 133 Main St., Afton, NY 13730 }
  153.